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Abstract. Wide variety of engineering design tasks can be 
formulated as constrained optimization problems where the 
shape and topology of the domain are optimized to reduce 
costs while satisfying certain constraints. Several mathemat¬ 
ical approaches were developed to address the problem of 
finding optimal design of an engineered structure. Recent 
works have demonstrated the feasibility of boundary 

element method as a tool for topological-shape optimization. 
However, it was noted that the approach has certain draw¬ 
backs, and in particular high computational cost of the itera¬ 
tive optimization process. In this short note we suggest ways 
to address critical limitations of boundary element method 
as a tool for topological-shape optimization. We validate our 
approaches by supplementing the existing complex variables 
boundary element code for elastostatic problems with robust 
tools for fast topological-shape optimization. The efficiency 
of the approach is illustrated with a numerical example. 

1 Introduction 

Problems of optimal design, i.e. variation of shape 
and topology of the domain to extremize certain functional 
subject to additional constraints, are ubiquitous in different 
branches of engineering. Recent emergence and rapid de¬ 
velopment of stereolithography and 3D printing technologies 
]3|4| have enabled fast prototyping of complex-shaped struc¬ 
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tures and revived the interest in automated optimal design. 
The most common formulation of optimization problems 
studied in stractural engineering (e.g. Q) seeks for an opti¬ 
mal shape and topology of an elastic body that minimize the 
strain energy (compliance) while satisfying the weight con¬ 
straint and the additional constraints imposed by the bound¬ 
ary value problem. Non-convex nature of such optimiza¬ 
tion problems often makes finding globally optimal designs 
difficult or impossible. Various numerical methods, includ¬ 
ing shape gradient-based approaches Q, level set methods 
homogenization 0[T§ and topological-shape sensi¬ 
tivity o have been developed during the last few decades 
to address this class of problems. These approaches are typ¬ 
ically implemented within the finite element method (FEM) 
context. In case if admissible designs include only homo¬ 
geneous regions with piecewise-constant properties (micro- 
structured composite designs are prohibited), the problem of 
optimization reduces to finding optimal configuration of the 
domain boundaries. Few recent papers 00 suggested that 
boundary integral approaches, and in particular boundary el¬ 
ement method (BEM), can be a convenient tool to address 
this class of problems. The implementation of optimization 
algorithms within BEM utilizes the concept of topological 
derivative (TD) | [TT| - [T3| - a cost of making an infinitesimal 
circular (spherical) cavity with a center in a given point of 
the domain. Existence of this kind of derivative opens wide 
avenue for numerous gradient-based approaches. The most 
straightforward one utilizes the idea of calculation of the 
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tion of topological derivative mi- Following p6| we define 
the TD as (Fig. 1(a)); 


D{x) = lim 
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Fig. 1. (a) Details on the definition of TD. (b) Optimization on 
the rectanguiar mesh, (c) Quadtree strategy of sampling of the 
TD inside the domain. 


field of TD on a given mesh, and removing material in the re¬ 
gions where the TD is below certain cutoff level. This kind of 
approach has demonstrated its feasibility in both 2D and 3D 
elastostatic problems Q|2). However, existing papers on the 
subject do not address the issue of complexity of topological- 
shape optimization with BEM, and, particularly, the question 
of competitiveness of BEM approaches as compared to EEM. 

In this note we suggest a simple optimization technique 
based on 2D BEM. It is demonstrated that even without using 
fast BEM techniques, one can easily reach 0{n^) asymptotic 
performance of shape optimization with BEM, which corre¬ 
sponds to the performance of 2D EEM approaches - (9(nv) 
{tis and Hy are the numbers of surface elements and volume 
elements after BEM/EEM discretization). 

The suggested technique is based on complex vari¬ 
ables boundary element method (CVBEM) It utilizes 
quadtree-based calculation of TD inside the domain, allow¬ 
ing 0{n^) performance of evaluation of helds inside the do¬ 
main. We also use simple 0{n^) algebraic solver based on 
blockwise update of the inverse system matrix to solve direct 
boundary value problem at every iteration. 

The paper has the following organization. The next sec¬ 
tion describes the method of topological optimization em¬ 
ployed in our research, as well as its novel features. Section 
3 provides an illustrative numerical example that is served to 
highlight the capabilities of our approach. The summary and 
discussion of our results, as well as the major directions of 
future work are presented in Section 4. 


2 Method 

2.1 Optimization technique 

In this paper we develop our approach in 2D using the 
framework of CVBEM | [T4l . The method offers a number of 
attractive features: i) complex hypersingular boundary inte¬ 
gral equation (BIE) for tractions and displacement discon¬ 
tinuities, derived for a system that can include an inhnite 
matrix, hnite-sized blocks with piecewise-constant proper¬ 
ties, voids and cracks ID; ii) closed form analytical calcu¬ 
lation of all the integrals in BIE | fT4| , iii) circular elements 
to approximate curved boundaries |14); vi) asymptotic ap¬ 
proximations to model cracks |[T4) ; v) implementation of 
piecewise-constant body forces |15| . 

Our scheme of topological optimization utilizes the no- 


Here 'F is the cost functional, fig is the original domain 
perturbed by the presence of an inhnitesimal cavity of ra¬ 
dius e centered at the point x, 5e is the small perturbation of 
the cavity radius, and / is regularizing, problem-dependent 
function. Here and below we’ll be working with the strain 
energy functional, that serves as a global measure of compli¬ 
ance of an elastic structure. Eor this case the area (volume) 
of the cavity serves as regularizing function. In absence of 
body forces, the closed form analytical expressions are avail¬ 
able both in 2D and 3D ||5][T^ . In the case of 2D plane strain 
elasticity the expression for TD writes 
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Here a and e are stress and strain tensors at point x. In 
case of presence of a constant body force b the expression (2) 
should be enriched with an additional term, associated with 
work done by the body force —bu{x), where u{x) is displace¬ 
ment at point X. The treatment of constant body force within 
CVBEM is discussed in 1151. It worth noting here that simi¬ 
lar analytical expressions for TDs are also available for other 
important classes of cost functionals, in particular, for the 
components of homogenized elasticity tensor of a periodic 
cell ifTTlITSl. 


Shape optimization approaches based on TDs use either 
hard-kill methods ED. or bidirectional optimization ID- 
In our work we utilize a hard-kill algorithm, in which the 
optimization is achieved by progressive elimination of mate¬ 
rial in the areas where the TD is below the cutoff level, as 
discussed further. 


2.2 Boundary generation procedure 

The new set of domain boundaries is created at every 
iteration of optimization procedure. Boundary generation 
strategy employed in our work is similar to one described in 
Q. We discretize the initial domain onto a set of M square 
cells (Eig. 1(b)). The helds inside the domain are calcu¬ 
lated at the center of each cell Xc{i), i = 1..M. The cutoff 
level for the TD is dehned as Dq — Dmin + Ot{Dmax — Dmin), 
where Dmin and D,nax are minimum and maximum values of 
TD within the current domain, and the coefficient a is tuned 
in range between 0.1% and 2%. Eor every cell of the ini¬ 
tial domain we dehne the Boolean status s. At the begin¬ 
ning of the iterative process s = 1 for every cell. At every 
iteration we assign s = 0 for the cells with D{xc{i)) < Dq 
(marked with empty points in Eig. 1(b)), and for isolated 
cells (marked with gray points in Eig. 1(b)). When the status 
































is assigned, we generate the new boundary using straightfor¬ 
ward algorithm - if i-th cell has i = 1 and its right (top, left, 
bottom) neighbor has s — 0 , then generate a corresponding 
boundary element. For every boundary between the neigh¬ 
boring cells one straight element with three points of collo¬ 
cation (quadratic approximation) is used. At every iteration 
we mark the elements that were deleted and those that were 
added at the current step. 

Boundary value problem constraints are incorporated 
explicitly - for the cells bounded by the elements with non- 
homogeneous Neumann boundary conditions s = 1. The vol¬ 
ume constraint does not present explicity in this scheme, so it 
is incorporated as the stopping criteria - the iterative process 
is discontinued when the ratio of current area of the material 
to the initial one R reaches the prescribed value Rq. 

2.3 Quadtree sampling of fields inside the domain 

The described procedure of sampling of TDs on the uni¬ 
form mesh should be considered computationally inefficient, 
since the calculation of a field of TDs inside the domain 
takes 0{n^) operations. In order to reduce the asymptotic 
complexity of this step, we use the specific strategy of cal¬ 
culation of TDs, which is based on quadtree algorithm of 
sampling p0| . The main idea is to use few different lev¬ 
els of refinement, that are employed for initial detection and 
further refinement of features of optimized domain. Within 
this approach each refined cell is subdivided onto 4 sub-cells 
(Fig. 1(c)). One can formulate different possible criteria 
of refinement. In this work we use the following criterion 
- if the values of s for a current cell and its nearest neigh¬ 
bor are not the same, both cells should be refined to the next 
level. The boundary generation algorithm described above 
is used to generate boundaries around the points of highest 
level of refinement. It is easy to see that such an algorithm 
of sampling and boundary generation requires calculation of 
the topological derivative at 0{ns) points inside the domain, 
leading to 0{n^) operations of evaluations of boundary inte¬ 
grals for each boundary element. 

2.4 Iterative update of the inverse matrix 

CVBEM | [T4) generates non-symmetric and non-sparce 
system matrix, and the corresponding system of equations 
is difficult to solve using regular iterative approaches. Here 
we describe one possible way to treat the resulting system 
matrix during the iterative updates of the initial boundary 
value problem. The original boundary consisting of el¬ 
ements leads to Ns = 6ns rows/columns in the resulting non- 
symmetric matrix generated by CVBEM (considering 3 col¬ 
location points per element and 2 degrees of freedom per col¬ 
location point). Assume that at k — th iteration elements 
have been added to a boundary, and elements have been 
removed. This results in Na = and W = added and re¬ 
moved rows/columns in the system matrix. If, as always the 
case in the iterative process, na ^ ns, n^ ^ ns, it is efficient 
to perform update of the inverse matrix using incremental 
expressions, such as Sherman-Morrison-Woodbury formula 
||2T1 or Banachiewicz ||22| formula for blockwise matrix in¬ 


verse: 


f ^ ‘ _ /A-' +A-^BS-^CA-^ -A-^BS-^ \ ,,, 

\Cd) “ -S-^CA-^ S-^ ) 

Where S = D — CA^^B. Denoting blocks of an extended ma¬ 
trix as E,F,G and H, and expressing A^\ we obtain the for¬ 
mula for inverse matrix update in case of removing elements: 


(4) 

It is clear that calculations of both expressions do not re¬ 
quire {Ns X Ns){Ns X Ns) matrix multiplications, only lower- 
rank operations with {Ns x W), {Ns x Na),{Ns x Nd),{Na x 
Na) and {Nd x Nd) matrices. This leads to 0{n^) perfor¬ 
mance of the algorithm of iterative matrix update. 

It is important to mention that in case if the matrix up¬ 
date is associated with the new cavity in the domain, matrix 
S is singular and requires regularization, e.g. truncated SVD 
regularization p3| (used in this work) or explicit incorpora¬ 
tion of additional conditions that fix rigid body motion of a 
new closed boundary G3- 

The suggested (or similar) schemes of low-rank inverse 
matrix update are well know and are widely used in adjacent 
fields of scientific and engineering calculations, including 
network structures, asymptotic analysis and boundary value 
problems with changing boundaries (see the excellent review 
presented in p4|). It is thus natural to adopt this approach for 
fast iterative optimization schemes with boundary elements. 

3 Numerical example 

In this section we discuss a simple benchmark example 
that allows us to evaluate the capabilities of our approach. 
Consider a 2D problem of optimization of a shape and topol¬ 
ogy of a fixed support. The initial domain has the square 
shape (Eig. 2(a)). The left side of the support is rigidly fixed, 
and the point load is applied at the upper-right angle. 

We utilize the optimization scheme described above to 
find an optimal shape and topology of the support, providing 
the smallest strain energy E for the volume fraction Rq = 0.5. 
The cutoff parameter a was set to 0.003. The TDs are sam¬ 
pled on a grid containing 20 x 20 cells. The solution was 
found in 29 iterations. The final shape is given in Eig. 2(b). 
The obtained solution is in good agreement with the one ob¬ 
tained with EEM and BEM in the earlier works MB- Fig- 
2(c) gives the evolution of the strain energy functional dur¬ 
ing the iterative process. Eig. 2(d) presents the comparison 
of times spent during the iteration with full LU factoriza¬ 
tion (LU), and the blockwise inverse matrix update, includ¬ 
ing adding rows/columns according to (3) {BU+) and remov¬ 
ing rows/columns according to (4) (BU-). The boost in per¬ 
formance clearly depends on rank of the update (number of 



operations, which is much faster than the corresponding 
operations performed within FEM (both take (9(nv) oper¬ 
ations). In addition, FEM techniques require good qual¬ 
ity domain discretization, generation of which takes at least 
0{njn{ny)) operations. 

These considerations motivate the development of three- 
dimensional fast BEM-based framework for topological- 
shape optimization of an elastic domain. Its implementation 
can be based on the principles similar to those presented in 
this work, and extended with additional features; smoother 
boundaries generation, combined shape and topology opti¬ 
mization iterations, etc. Clearly, the BEM-based techniques 
should become the most computationally efficient tools for 
topological-shape optimization. 


Fig. 2. Shape optimization of a fixed support, (a) Initial bound¬ 
ary value problem, (b) Final configuration reached after 29 itera¬ 
tions. (c) The ratio of current strain energy E to the initial strain 
energy Eq as a function of current volume fraction R. 


added/removed elements) and varies significantly from iter¬ 
ation to iteration. The total time spent on the full LU factor¬ 
ization is 41 s, whereas the total time spent on the blockwise 
matrix update is 15.5 s. Calculations were performed with 
regular Core i-5 laptop machine. Linear algebra operations 
were implemented within non-optimized BLAS/LAPACK 
framework I 


and win-32 gfortran compiler 1261. 

In order to illustrate the quadtree sampling of TDs in¬ 
side the domain, we consider the same example, but with 
the twice higher degree of grid rehnement. Two levels of 
grid refinement are used. The coarse level discretizes do¬ 
main onto 20 X 20 cells. As we could see above, this level 
is sufficient to detect the important features of the optimized 
domain. On the finer level we discretize the domain onto 
40 X 40 cells. This level is used to render finer features of 
the optimal design. Fig. 2(e) gives the sampling points dur¬ 
ing first iteration. Fig. 2(f) demonstrates corresponding con¬ 
figuration of domain boundaries. Using quadtree sampling 
has led to decrease of the number of points inside the do¬ 
main from 1600 to 676, and decrease of the corresponding 
computational time from 11.3 to 4.8 s. 


4 Discussion and conclusions 

In this work we suggested a set of tools for topological- 
shape optimization with boundary elements. As we could 
see on the example of 2D CVBEM method and a simple 
toolkit for topological optimization, one can reach the com¬ 
putational performance available with 2D FEM techniques. 

It is therefore clear that the usage of fast BEM tech¬ 
niques 127 281 would undoubtedly allow to outperform all 
the existing FEM techniques of topological-shape optimiza¬ 
tion. For example, using fast BEM would allow solution of 
direct boundary value problem for 0{ns) operations, and cal¬ 
culation of the field of TD inside the domain for 0{njn{ns)) 
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